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ABSTRACT 

Recently,  Micchelli,  Rivlin  and  Winograd  [ 15]  described  a scheme 
for  interpolating  by  splines  at  given  data  points  which,  in  a certain  reason- 
able sense,  is  optimal  among  all  schemes  which  attempt  to  recover  a func- 
tion from  its  values  at  those  data  points.  This  paper  offers  a Fortran  program 
for  the  calculation  of  that  optimal  interpolant.  A short  derivation  of  that 
recovery  scheme  is  given  first,  in  order  to  make  the  paper  selfcontained  and 
also  to  provide  an  alternative  to  the  original  derivation  in  [ 15].  For  the 
latter  reason,  a derivation  of  the  related  envelope  construction  of  Gaffney 
and  Powell  [ 8]  is  also  given.  From  a computational  point  of  view,  these 
schemes  are  special  cases  of  the  following  computational  problem:  to 
construct  an  extension  with  prescribed  norm  of  a linear  functional  on  some 
finite  dimensional  linear  subspace  to  all  of  IL  [a,  b]. 
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1.  INTRODUCTION 

This  paper  offers  a Fortran  program  for  the  calculation  of 
the  optimal  recovery  scheme  of  Micchelli,  Rivlin  and  Winograd  [15]. 
A short  derivation  of  that  recovery  scheme  is  given  first,  in 
order  to  make  the  paper  self contained  and  also  to  provide  an 
alternative  to  the  original  derivation  in  [15].  For  the  latter 
reason,  a derivation  of  the  related  envelope  construction  of 
Gaffney  and  Powell  [8]  is  also  given.  From  a computational  point 
of  view,  these  schemes  are  special  cases  of  the  following  computa- 
tional problem:  to  construct  an  extension  with  prescribed  norm 

of  a linear  functional  on  some  finite  dimensional  linear  subspace 
to  all  of  ILj[a,b]. 


2.  THE  OPTIMAL  RECOVERY  SCHEME  OF  MICCHELLI,  RIVLIN  AND  WINOGRAD 


This  scheme  concerns  the  recovery  of  functions  from  partial 


information  about  them  Let  n > k and  let  t 


decreasing,  in  some  interval  [a,b],  with  t.  < t . , 

. . _ (k)  , „(k-l) . ..  (k-1) 

For  f f l'  :»{geC  [a,b]  : g 


be  non- 

all  i. 
(k) 


( IL  },  denote  by  f 
00 

sequence  x,  i.e.,  ) 


the  restriction  of 


abs.  cont . , g 
f to  the  data  point 


:>  (f  )n 

vVi 


with 


J j (i)  max(m  : 


Ti} 
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! (k)  (k) 

We  call  a man  S : IL  *►  IL  a recovery  .scheme  (on  I.  ) 

I QO  OO 1 OO 

with  respect  to  t provided  Sf  depends  only  on  fj^, 
provided  f j ^ = gj^  Implies  that  Sf  = Sg.  The  (possibly 

infinite)  constant 

constg  :=  s up { ||  f — Sf  ||  / ||  f ^ ||  : f c n,  ^ } 

measures  the  extent  to  which  such  a recovery  scheme  S may  fail 
to  recover  some  f since  it  provides  the  sharp  error  bound 

(1)  ||  f - Sf  ||  < const,,  ||  f (k)  ||  , all  f e 1L  <k)  . 

Here  and  below. 


||  g ||  :=  ess.  sup  { |g(t)|  : a <.  t ^ b} 


A recovery  scheme  S is  optimal  if  .its  const  is  as  small 

a 

as  possible,  i.e.,  if  the  worst  possible  error  is  as  small  as 
possible,  as  measured  by  (1).  We  write 

const  :*  inf  const 
T S b 

for  the  best  possible  constant.  Here  is  a quick  lower  bound  for 
that  constant.  We  have 


const_  “ sup  sup 

f g 


g ~ Sf  ||  /||  g 


(k) 


sup  ||  g - S0||  /||  g 


(k) 


- sup  max{  |i  g - SO  || , ||  -g  - SO  ||}/  ||  g^  || 

e|t-° 

1 sup  II  g||/||  g<k>  II  , 

g|x“° 

for  every  recovery  scheme  S , hence 


II  sll/ll  8<k>  II  - 


(2)  const  > c(t)  :=  sup 

sl*c 


-1- 


I 


I 


But,  actually,  equality  holds  here.  For  the  proof,  we  need  the 
notion  of  perfect  splines. 

A perfect  spline  s of  decree  k with  (simple)  knots 
< . . . < 5 in  [a,b]  is  any  function  s of  the  form 


s(x)  = P(x)  + c 


r , .1  r i+1,  .k-1  . 

I (~)  / (x  - y)+  dy  (c 


i=0 


0"' 


^r+l ‘=  b) 


with  P e : = polynomials  of  decree  < k and  c some  constant. 

In  other  words,  such  a function  s is  any  k-th  anti-derivative 
of  an  absolutely  constant  function  with  r sign  changes  or  iumps 
in  [a,b].  Their  connection  with  the  present  topic  stems  from 

Pro£os^t_ion  1 (Karlin  [10]).  s is  a perfect  spline  of 

degree  k with  n - k knots  then 

||s(k,||  ■ inf(||f(t)||  : f . n .<k>,  f|t  . .|T>  . 

This  proposition  follows  directly  from  the  following  lemma. 

(k) 

Lemma  1.  If  h c L is  such  that,  for  a = < ...  < C 

■■■  — °°  U r+1 

* b. 

(k) 

a^h  > 0 a . e . on  for  some  c {-1,1},  i=0,...,r  , 

and  h | ^ ■ 0 , then  r ^ n - k.  Further,  if  r n - k,  hence 
r * n - k,  then 

(3)  Ti  " Ci  < Ti+k’  ^ 1 • 


Assuming  this  lemma  for  the  moment,  we  see  that  if,  in 

Proposition  1,  we  had  |j  s^||  > ||  f^^||  for  some  f € IL^^  with 

f|T  * s | t » then  h :*  s - f would  satisfy  the  hypotheses  of  the 

lemma  for  some  r < n - k,  which  would  contradict  the  conclusion 
of  the  lemma. 

Proof^fiiiLemma  1:  By  Rolle's  theorem,  there  must  be  points 


C1  < 


< t 


n+l-k 


(k-1) 

at  which  h vanishes, 


and  for  which 


y 


We  are  now  fully  equipped  for  the  attack  on  the  optimal 
recovery  scheme.  By  Karlin's  theorem,  the  set 

Q :*=  {q  : q is  a perfect  spline  of  degree  k with  £ n - k 
knots,  q I =0} 

is  not  empty.  By  Lemma  1,  each  q e Q^\{0}  has,  in  fact,  exactly 

n - k knots  and  does  not  vanish  off  t (since  otherwise  it  would 
have  n + 1 zeros  yet  only  n - k sign  changes  in  its  k-th 

(k) 

derivative,  contrary  to  the  lemma).  Therefore,  if  f e L^  and 

f|T  ■ 0 and  x / t, then  (f(x)/q(x))q  is  a well  defined  perfect 

spline  of  degree  k with  n - k knots  which  agrees  with  f at 
the  n points  of  t and  the  additional  point  x,  hence, 
Proposition  1 implies  that 

(5)  ||f(k)||  > | f (x)/q(x)  | ||  q°°||  . 

It  follows  that 


(6)  sup(  | f (x)  | / ||f (k>  !|  : f e l/k),  f|T  = 0}  ■=  |q(x)|/||  q(k)||  . 

Since  the  left  side  is  independent  of  q,  this  shows  that  is 

the  span  of  just  one  function,  say  of  q,  normalized  to  have 

q(k)(0+)  = 1 . 

It  also  shows  that 

(7)  c(t)  = ||  q || 


This  gives  a means  of  computing  c(t),  but  not  quite  yet  the 
equality  const  = c(t)  nor  the  optimal  recovery  scheme.  For  this, 

let  £,,...,£  be  the  n - k knots  of  q and  consider 
1 n-k 

$ :=  splines  of  order  k with  simple  knots  £ , . . . ,£n_k 

- ^(k-2) 


:=  IP,  _ n C 


: = { f e C 


(k— 2 ) 


- . e P , all  i}  U0-a,£ 

’i’^i+1' 


n+l-k 


Theorem  2 (Micchelli,  Rivlin  & Winograd  [15]).  The  rule 

(ff)  = f|  and  Sf  e $,  all  f e IL  (k) 

A (k) 

defines  a map  S (read  "ess  crown(ed)")  on  IL  which  is 

linear  and  an  optimal  recovery  scheme  with  respect  to  x . 

Proof.  First,  we  prove  that  ^ is  well  defined.  Since,  by 
the  lemma,  x < £ < ^i+k>  we  ask 

for  matching  f ^k  ^ at  some  point  x only  when  such  x / £, 

/i  -I  \ 1 

i.e.,  when  s^  ;(x^)  makes  sense  for  each  s e Secondly, 

dim  $ = k + number  of  polynomial  pieces  - 1 = n,  therefore  S is 
well  defined  if  and  only  if 

(8)  f j ^ =■=  0 and  f e $ implies  f = 0 . 

For  this,  we  would  like  to  use  inequality  (5),  but,  since  $ is 

(k) 

not  contained  in  IL  but  only  in  the  larger  space 


=-  " nn  «-[k)«i-W  • 

i«0 


I 


we  must  first  extend  (5)  to  such  spaces,  with  the  convention  that 


00 


max 

i 


■ 00 


^i’^i+l5 


for  fell. 


00 


This  requires  the  following  slight  strengthening  of  Lemma  1. 

(k) 

I.emma  1'.  JLf  h e IL  ^ for  some  a = < . . . < £ 

and, for  some  o e {-1,1},  r+ 


o(-)xh 


i.00 


> 0 a. e.  on  (C  ± , C » 1 = 0,...,r  , 


then  h | ^ =0  implies  n - k r . 

Here,  we  mean  by  h^  ^(x)  = 0 that  h^  ^ (x  )h^  ^(x+)  <_  0, 
in  case  the  number  x occurs  k-fold  in  x. 


Proof.  Rolle's  theorem  gives  again  n + 1 - k distinct 
(k-1) 

zeros  of  h - which  again  consists  of  r + 1 strictly  monotone 

pieces,  but  may  fail  to  be  continuous  across  the  points  £ . 

(k)  1 r 

But,  since  h alternates  in  sign,  this  latter  fault  can  easily 

be  remedied  by  appropriate  local  linear  interpolation  across 


a small  neighborhood  of  each  discontinuity  without  disturbing 

the  other  two  properties  and  now  n - k r follows  as  before;  Q.E.D. 

This  gives 

Proposition  1'.  If  s is  a perfect  spline  with  < n - k 
knots,  say  with  knots  < ...<  t where  r < n - k,  and  of 


degree  k,  then  even 

,(k)  ||  _ j„ r / II  r (k) 


l|s'K'||-  intill  r"'ll  : f « U-i1;,  Cl 

From  this  we  conclude,  with  ||q^^||=  1.  that 
(5')  ||  f(k)||  > |f (x)/q(x)  | for  x i T,  f € l/k>,  f|T  = 0 


s } . 

T 


But  this  implies  (8)  since 


(k) 


0 for  f c $,  and  so  shows 


that  S is  well  defined.  Finally,  (5')  also  implies  that,  for 

f n <k> 
f e ^ » 

||  f (k)  ||  - ||  (f  - Sf)(k)||  > |f(x)  - £f(x)|/|q(x)| 


or 

(9)  | f (x)  - Sf(x)J  < |q(x)|  ||  f(k)|| 

showing,  with  (2)  and  (7),  that  f is  an  optimal  recovery  scheme. 
This  proves  Theorem  2. 

(k) 

MRW  actually  insist  that  a recovery  scheme  S map  into  IL^  , 
hence  they  are  not  quite  done  at  this  point,  since  S only  maps 

into  lL^k  But,  since  $ can  be  viewed  as  spline  functions 

of  degree  k with  double  knots  at  the  £ , we  can  produce  an 
(k) 

element  of  IL^  arbitrarily  close  to  Sf  merely  by  pulling  all 

these  double  knots  apart  ever  so  slightly.  This  shows  that 
inf  const.  c c(t)  even  if  the  inf  is  restricted  to  S mapping  into 
(k)  s 

IL  but  now  the  inf  is  not  attained  apparently  for  k >1. 


3.  THE  ENVELOPE  CONSTRUCTION 

The  preceding  discussion  allows  a simple  derivation  of  sharp 

(k) 

estimates  for  the  value  of  a function  f e IL  at  some  point  x, 

00 

given  the  vector  f|  and  a bound  L on  its  k-th  derivative  on 
[a,b],  as  follows.  ' 

We  are  to  construct  the  set 

I (f(x)  : f e F} 
x 

with 


7 


F {f  £ 1L  <k)  : f|T  - a,  ||  f(k)||  < L} 

for  some  given  a and  L.  If  F is  not  empty,  then  I is  a 
closed  interval, 

I [a  ,b  ] 

X 1 X*  X 

say,  since  F is  closed,  convex  and  bounded  and  [x]  : f w-  f (x) 

00 

is  a continuous  linear  functional  on  1L  . Assume  that 

oo 

F°  : “ { f c F : ||  f(k)  1 1 < L } * 0 . 

Then 

(10)  (ax,bx)  - [x]F°. 

Karlin's  theorem  then  implies  that,  for  x / x, 

:»  {q  : q is  perfect  spline  of  degree  k with  £n  - k 
knots,  q|T  - a,  q(x)  = ax) 

is  not  empty.  Further,  Qx  £ F,  since,  by  definition  of  ax> 
there  exists  f £ F with  f(x)  = ax  while  each  q e agrees 
with  such  an  f at  the  n + k + 1 points  x and  x,  therefore 
II  q^||£  ||  f^||  £ L,  by  the  proposition.  On  the  other  hand, 

0 n F°  - 0,  since,  if  q e 0 n F°,  then  a ■ q(x)  £ (a  ,b  ) , 

^ X X XX 

by  (10),  a contradiction. 

It  follows  that 

for  all  y,  ay  £ q(y)  £ by  . 

But,  for  any  g£F°,  h :■  g - q has  £ n - k sign  changes  in  its 
k-th  derivative  and  vanishes  at  x,  therefore  q has  exactly 
n - k knots  and  h does  not  vanish  off  x,  by  Lemma  1.  Hence, 
if  Ti+i» * * * ,Ti+j  are  t*ie  points  of  x between  x and  y i x, 
then 

(-)^(g(y)  - q(y))  > o . 


This  shows  that 


q(y) 


if  the  interval  between  x and  y contains  an 


-8- 


number  of  points  of  x.  It  follows  that  Q contains  exactly  one 

function  and  that  this  function  supplies  half  of  the  entire 
boundary  of 

(11)  {(x,f(x))  : x e [a,b],  f e F}  , 


the  other  half  being  supplied  by  the  perfect  spline  p of  degree 

k with  n-k  knots  for  which  p|  = a and  p(x)  = b . 

x x 


We  note  the  curious  corollary  that  the  perfect  spline  s of 
Karlin's  theorem  is  unique  if  there  exists  g which  agrees  with 

f at  all  points  of  x but  one  and  for  which  ||  ||  < ||  s^^||  , i.e. 

if  at  least  one  of  the  interpolation  constraints  is  active.  Put 
differently,  it  says  that  if  there  are  two  different  perfect 
splines  s,  s of  degree  k with  _<  n - k knots  for  which 
si  = si  = a and  which  agree  at  some  x / x,  then 


l|s(k)||-  II 


;(k,ll- 


La  :=  min{  |i  f 


(k) 


f £ L 


(k) 


a} 


Let  now  q be  the  half  of  the  envelope  of  (11)  with 

(v)  4-  fit')  4- 

q 7 (0  ) = L,  and  let  p be  the  one  with  p ' (0  ) = -L.  Gaffney 
and  Powell  [8]  choose 


SLa  :=  (p  + q ) / 2 

as  a good  interpolant,  its  graph  being  clearly  the  center  of  (11). 
Since  p and  q are  uniquely  defined  for  L > by  the  requirement 

that  they  are  perfect  splines  of  degree  k with  n-k  knots,  equal 
to  a at  x and  t.o  a^,  resp.  b^  at  x,  they  are  necessarily  continuous 

functions  of  L and  a in  that  range.  In  particular,  with  q =:q  , 

L y 01 

we  have  qT  = Lq,  and  q,  , + q,  . = o as  L + »>,  and, 

L,a  ^l,a/L  l,a/L  ^1,0 

similarly,  p , + -q  as  L + ®.  This  shows  that,  with 

L>  f 01  / La 

u,  < . . . < u . the  knots  of  q and  v,  < . . . < v , the  knots  of  p, 
1 n-k  n 1 n-k 

lim  u ■ lim  v = £ , i = l,...,n  - k . 

L-h»  1 L-*»  1 1 


In  particular,  for  large  enough  L,  the  sequence 
0 = C0+,C1-,C1+* ‘ ' ' ,Cn-k-,Cn-k+’^n-k+l-  “ 1 

with 

:■  min{ u^ , v^ } , Ci+  t*  max{u^,v^} 


is  nondecreasing  and 


<SI.a) 


(k) 


f ° °n 

\±2L  on  (C±_,Ci+) 


hence 

II  <4 SL«)  (k>  II  i 1 2L||  u - v|| ► 0 . 

L-*°° 


This  shows  that  S^a  converges  to  an  element  of  $,  while 

S. a | = a for  all  L , 

L [ T 

therefore,  S^a  converges  to  the  optimal  interpolant  for  the  data. 

It  was  in  this  way  that  Gaffney  and  Powell  [8]  constructed, 
quite  independently  from  Micchelli,  Rivlin  and  Winograd  [15],  the 
same  optimal  recovery  or  interpolation  scheme  S. 

The  problem  of  constructing  the  set  I was  posed  originally 

in  the  basic  paper  by  Golomb  and  Weinberger  [9],  although  they 
gave  detailed  attention  to  such  problems  only  when  the  (semi)norm 
involved  comes  from  an  inner  product.  Micchelli  and  Miranker  [14] 
solved  the  problem  posed  at  the  beginning  of  this  section  in  the 
sense  that  they  correctly  described  the  boundary  of  (11)  as  being 
given  by  just  two  perfect  splines  of  degree  k,  each  with  n - k 
knots,  and  with  their  k-th  derivative  equal  to  L in  absolute  value. 
In  fact,  Micchelli  and  Miranker  consider  the  slightly  more  general 
(k) 

situation  where  f is  only  known  to  map  [a,b]  into  some  interval 
[m,M],  They  state  that  these  matters  could  be  proved  along  the 
lines  used  by  Burchard  [6]  to  solve  a related  restricted  moment 
problem  and  refer  specifically  to  Karlin  and  Studden  [11,  VIII,  53] 
for  requisite  facts  concerning  principal  representations  of  interior 
points  of  moment  spaces.  Of  course,  these  facts  go  back  to 
Krein  [12].  Quite  independently,  Gaffney  and  Powell  [8]  also 
solved  this  problem,  with  the  proofs  in  Gaffney's  thesis  based  on 
Chebyshev  type  inequalities  as  found  in  Karlin  and  Studden  [11, 

VIII,  §8]  and  adapted  by  him  to  weak  Chebyshev  systems. 


4.  THE  CONSTRUCTION  OF  NORM  PRESERVING  EXTENSIONS  TO  ALL  OF  IL 

A 

Both  the  optimal  recovery  scheme  S of  Section  2 and  the 
envelope  of  Section  3 require  the  construction  of  an  absolutely 
constant  function  h with  no  more  than  a specified  number  of  jumps 


-10- 


which  provides  an  integral  representation  of  a linear  functional 
given  on  a linear  space  of  spline  functions. 


In  the  optimal  recovery,  we  are  to  construct  a perfect  spline 

s of  degree  k with  n - k knots  and  with  ||  s^k^||  = 1 which  vanishes 

at  the  given  n-point  sequence  t.  Let  (M^)"  be  the  sequence  of 

B-splines  of  order  k with  knot  sequence  t,  each  normalized  to  have 
unit  integral,  i.e.. 


(12,  Hl(x)  Mlk>t 

with  hit . . . .Ti+klf 
f at  the  points  t^,. 
integral  remainder, 


(x)  :=  k [ t i , . . . , x i+k ] (*  - x)k  1 , 

the  k-th  divided  difference  of  the  function 
’*,Ti+k'  ^en»  from  Taylor's  expansion  with 


It 


i... 


i+k 


Jf  - / i+kM  . (x)f(k)(x)dx/k!  for  all  f 

1 , X , T 


e I 


(k) 


n k 

The  points  £ = are  therefore  characterized  by  the  require- 

ment that  the  function 


n-k 


h (x)  :=  sign  "]  [ (x  - E) 

S j _ 1 *■ 


= 


(x) 


i=l 


be  orthogonal  to  each  of  the  n-k  functions 


Before  considering  the  computational  details  of  determining 
£ from  this  orthogonality  condition,  I want  to  comment  on  the 
fact  that  this  is  a problem  of  representing  or  extending  a linear 
functional  on  some  subspace  of  IL  ^ and  is  therefore  closely  related 

to  the  problem  of  computing  the  norm  of  a linear  functional  on  some 
subspace  of  IL  . This  is  also  explored  in  a forthcoming  paper  by 
Micchelli  [13]. 

Indeed,  if  T is  a linear  subspace  of  IL^  of  dimension 
n + 1 - k,  and  X is  a linear  functional  on  T,  we  might  ask  for  £ 
and  a so  that 

(13)  a / h^g  = Xg  for  all  g c T . 

But  then,  in  particular,  h^  is  orthogonal  to  ker  X,  a subspace  of 

dimension  n-k.  Conversely,  If  we  have  already  found  h^  ortho- 
gonal to  ker  X then  there  will  be  exactly  one  a so  that  ah^ 
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represents  X on  T in  the  sense  of  (13),  unless  is  even  ortho- 
gonal to  all  of  T.  But  this  latter  event  cannot  happen  in  case  T 
is  weak  Chebyshev  since  h^  has  only  n - k jumps. 

It  is  clear  that  any  such  representation  h^  for  X produces 

an  upper  bound  for  the  norm  of  X.  In  fact,  ||  X||  ■ |a|  in  case  T 
is  weak  Chebyshev.  This  is  actually  how  I became  Interested  tiro 
years  ago  in  the  numerical  construction  of  representers  of  linear 
functionals  [3],  [4].  I was  interested  in  computing,  or  at  least 
closely  estimating,  the  norm  of  certain  linear  functionals  on 
certain  spline  subspaces  in  IL^.  In  a way,  this  is  a trivial 

problem,  viz.  the  maximization  of  a linear  function  over  a finite 
dimensional  compact  convex  set,  and  there  was  the  feeling  that 
there  ought  to  be  special  methods  available.  Perhaps  some  reader 
can  steer  me  towards  such  methods.  I found,  for  the  particular 
cases  of  concern  to  me  in  which  T was  always  weak  Chebyshev, 
nothing  more  effective  for  calculating  ||  x||  than  to  construct 
such  a representation  (13). 

Finally,  the  envelope  construction  corresponds  to  the  slightly 
different  situation  where  dim  T * n - k,  X e T'  and  a with 
| a | > ||  X ||  is  prescribed  and  one  seeks  £ so  that  again  (13)  holds. 
We  have  now  one  less  condition  to  satisfy  but  also  one  less  para- 
meter to  do  it  with. 


5.  CONSTRUCTION  OF  THE  KNOTS  FOR  THE  OPTIMAL  RECOVERY  SCHEME 

n_k 

As  we  saw  in  the  preceding  section,  the  knots  £ * (£^)" 

for  the  optimal  recovery  scheme  are  the  solution  to  the  following 

problem.  We  are  given  (Ti)"  nondecreasing,  with  < T^+jt»  all  i, 

and  n > k,  and  wish  to  construct  £_<...<£  with 
* 0 *r 

r : - n - k + 1 

and  with  £_  » a :■  x, , £ ■ b :■  t so  that 

u x r n 

T £ . 

(14)  l (S.  f 3 M.  - 0,  i - 1 n - k 

j-1  J ^ 

while  also 


(15)  Bj  + BJ+1  - 0,  J - 1 n - k . 
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Extend  t by 

n n+1  n+2  n+k 

Then  one  verifies  easily  that 


/ Mi,k  = Nm, k+1 (x)  f°r  X 


m>i 


with 


N 


Tm+k+l  Tm 


tJ.t.  a11  " • 


m,k+l  ' k+1  m,k+l 

Therefore,  (14)  is  equivalent  to 

jl  - N..W(t)-l,)  ’ °-  1 ■ l--"  - k • 

or,  on  subtracting  equation  i from  equation  i - 1 for  i = 2,...,n-k, 

W+l'V  - Bi,k+l<Vl»  • °’  1 ■ 1 " - k - 1 

l l OVm'V  - Nm,k+l<£j-l»  • 0 • 

j=l  J m>n-k  * j » j 

Using  the  fact  that  N ,,,(£_)  = 0 for  m > 1,  this  can  also  be 
. m, K > 1 0 

written 

(Bj  ~ 6j+l)Ni,k+l(^j)  = ~6rNi,k+l(Cr)’  1 = “ k " 1 

n-k 

l ~ 6,,.)  I N - -6  l N (E  ) • 

.L.  i 1+1  ^ , m,k+l  j r , m,k+l  r 

j-1  J J m>n-k  J m>_n-k 

Since  k+j_(£r)  = 0 for  i < n - k,  while  £ k+1(Cr)  * 1, 

the  right  side  becomes  simply  (0,...,0,-8.  ?>n-k 

Choose  now  8 = -1  to  make  things  definite.  Then  8 *=  (~)r  * 

n_t_ 4 r J 

■ (-)  J by  (15),  and  (19)  and  (15)  are  seen  to  be  equivalent  to 

(16)  F(U  - 0 

with  F : lRn*c-*Rn*c  given  by 


••  I >- 


i < n - k 


(17)  F(C)i  < 


5"  a , 1 * n - k 

“ m 


m=n-k 


where 

nyk  (_)n-k-j 

dS)  a.  ^ ji,  ( > 

-1/2 


N1 ,k+l ^ ’ 1 * " 1 

, 1 «■  n . 


We  solve  (16)  by  Newton's  method.  From  the  current  guess  £,  we 
compute  a new  guess  £*  = £ + 5£,  with  the  solution  to  the 

linear  system 

(19)  F'(06C  = -F (£)  . 

Since  N|  = Mi  k ~ ^i+1  k’  at^ition  equation  i to  equation 
i-1,  i = n - k,...,2,  in  (19)  produces  the  equivalent  linear  system 


1 l k-M  k)(c.)  (-)n_k"j5C.  = - l a , i=l, . . . ,n-k. 

j=l  m=i  m’k  m+1’k  ^ 3 m=i  m 

n-1 

But,  since  £ (M  - M ,)(t)  = M (t)  - M (t)  = M (t) 

£ m,k  m+l,k  i,k  n,k  i,k 

for  t <_  b,  this  shows  that  (19)  is  equivalent  to  the  linear  system 

(20)  Cx  = d 
with 

(21)  . v ■ »i>k<V 

i,j  = l,...,n  - k. 

The  matrix  C is  totally  positive  and  (2k  - l)-banded,  hence  can  be 
stored  in  2k  - 1 bands  of  length  n - k each,  and  can  be  factored 
cheaply  and  reliably  within  these  bands  by  Causs  elimination  with- 
out pivoting  (see  [5]). 

The  iteration  is  carried  out  in  the  program  SPLOPT  below, 
starting  with  the  initial  guess 

(22)  £<  - (t1+1  + ...  + t<xtl)/(k  - 1),  i ■ 1, . . . ,n  - k . 


i+k-1 


- 


A first  version  of  the  program  was  equipped  to  carry  out  Modified 
Newton  iteration:  £*  is  computed  as  the  first  vector  in  the 
sequence 


£ + 2-h6C,  h = 0,1,2, .. . 

for  which  ||  F(£  + 2 ^5£)  ||  ^ K II  F(C)||  2*  » 1°  all  examples 

tried,  the  initial  guess  (22)  was  apparently  sufficiently  close  to 
the  solution  to  have  always  h <=  0,  i.e.,  ||F(£  + 6£)  ||  2 < 11^(0  ||  2 
In  fact,  the  termination  criterion 

II  «ell  . 1 10"6(Tn  - V/(n  ' k) 


was  usually  reached  in  three  or  four  clearly  quadratieally  converg- 
ing iterations.  For  this  reason,  the  program  SPLOPT  below  carries 
out  simple  Newton  iteration.  It  would  be  nice  to  prove  that  Newton 
iteration,  starting  from  (22),  necessarily  converges.  But,  such  a 
proof  would  necessarily  have  to  be  in  control  of  the  norm  of  the 
inverse  of  the  matrix  F'(£),  hence  in  control  of  the  norm  of  the 
inverse. of  C = (N (£.))  as  a function  of  £.  Good  estimates  for 
— 1 ^ 

||  (N^(£^))  ||  have  been  searched  for  in  the  past  by  some  who  were 


interested  in  bounding  (the  error  in)  spline  interpolation,  but 
without  much  success.  E.g.,  the  simple  conjecture  that,  for  the 


initial  guess  (22),  ||  (N  (£.))  *||  <_  const.  Independent  of  £ has 

t J i K 

been  proved  so  far  only  for  k <_  4. 


SUBROUTINE  SPLOPT  < TAU.  N.  K»  SCRTCH.  T.  IFLAG  ) 

COMPUTES  THE  KNOTS  T FOR  THE  OPTIMAL  RECOVERY  SCHEME  OF  ORDER  K 
C FOR  BATA  AT  TAU  t I » . i*i . , . . » N . TAU  MUST  BE  STRICTLY  INCREASING. 
C SEE  TEXT  FOR  DEFINITION  OF  VARIABLES  AND  METHOD  USED. 

C IFLAG  = I OR  2 DEFENDING  ON  WHETHER  OR  NOT  T WAS  CONSTRUCTED. 

C DIMENSION  SCRTCHC < N-K > * ( 2*K  + 3 > +5  UK + 3 > > T(NIK) 

DIMENSION  TAU(N) . SCRT CH < 1 ) ,T< 1 ) 

DATA  NEWTMX.TOLRTE  / 10. . 000001/ 

NMK  = N-K 

IF  ( NMK ) 1.54.2 

1 PRINT  40 l.N.K 

601  FORMAT  < 13H  ARGUMENT  N = »M,29H  IN  SPLOPT  IS  LESS  THAN  K = ,I3) 

GO  TO  999 

2 IF  <K  .GT.  2)  GO  TO  3 

PRINT  402. K 

602  rORMAT < 13H  ARGUMENT  K =.I3»27H  IN  SPLOPT  IS  LESS  THAN  3) 

GO  TO  999 

3 NMKM1  * NMK  - 1 
FLOATK  » K 
KPK  « K+K 
KP1  * K+l 
KPKP1  « K+KP1 
KM1  = K-l 
KPKM1  « K+KM1 
KPN  * K+N 
SIGNST  * -1. 

IF  <NMK  .GT.  ( NMK/2  )*2 ) SIGNST  «=  1. 


-1 


fN4K4K 


C SCRTCM(l)  = TAH-rxTENDED(I>F  1-1,... 

NX  «-  N 4KFKP 1 

C SCRTCIKI4NX)  - Xl(l>.  1=0,  . . . ,N-mi 
NA  =•  NX  * NMK  * 1 

C SCRT CH(HMrt)  - -A(I)f  1-1,..., N 
NP  = NA  + N 

c scRTrn(HNti)  = xd)  or  dcd.  n-k 

NV  --  ND  4 NMK 

C SCRTCH( H MV ) = VNIKX(I)fI=1f...fK41 
NC  = NV  4 KF1 

r SCRTHU ( J-1)*(N-K;4I  4 NC)  = CHAT (I , J> f 1= 1 , . . . , N-K f J=1 f . . . . 2*K-t 
LENGCH  - NMKtKRKMl 

C EXTENT'  TAU  TO  A KNOT  SEUUENCE  AND  STORE  IN  SCRTCH. 

DO  5 J- 1 * N 

SCRTCH  ( J ) =-  TALK!) 

5 SCRTCH (KRN4J)  = TAU(N) 

DO  A J=1 ,N 

6 SCRTCH  ( K 4 .1 ) = TAU<J> 

C FIRST  GUESS  FOR  SCRTCH  ( • 4NX ) = XI  . 

SCRTCH  < NX ) - TAU<1) 

SCRTCH ( NMK414NX)  = TAU(N) 

DO  10  J=1»NMK 
SUM  = 0. 

DO  9 L = 1 f KM 1 

9 SUM  = SUM  4 TALKJ4L) 

10  SCRTCH ( J4NX ) = SUM/KM  t 

C LAST  ENTRY  OF  SCRTCH  ( . 4NA ) = - A IS  ALUAYS  ... 

SCRTCH  ( N 4NA  ) =-  .5 
C START  NT U TON  ITERATION. 

NEWTON  = 1 

TOT.  = TOLRTE*(TAU(N)  - TAU(1))/NMK 
C START  NEWTON  STFP 

COMPUTE  THE  2K-1  RANDS  OF  THE  MATRIX  C AND  STORE  IN  SCRTCH (. i NC > , 

C AND  COMPUTE  THE  VECTOR  SCRTCH (. 4NA ) - -A. 

00  DO  r*S  1 = 1 ,L£HGC.H 

21  SCRTCH ( 1 4NC  > = 0. 

DO  2?  I =2. N 

22  SCRTCH  < I - 1 4NA ) = 0. 

SIGN  = SIGNST 

ILEFT  = KF1 
DO  28  J= 1 . NMK 

XIJ  = SCRTCH ( J4MX  > 

23  IF  <XIJ  .LT.  SCRTCHT ILEFT41  ) ) GO  TO  25 
ILEFT  = HEFT  4 1 

IF  (ILEFT  .LT.  KFN)  GO  TO  23 

ILEFT  *•  ILFFT  - 1 

25  CALL  PSRL  VN ( SCRTCH  fK.IfXIJf ILEFT  »SCRTCH(1 4NV  ) ) 

ID  = MAXO (0. ILEFT-KPK ) 

INDEX  * NC4< J-ID4KM1 >*NMN+ID 
LLMAX  = MINO(KfNMK-ID) 

LLMIN  = 1 - MINOLOf ILEFT-KPK) 

DO  26  LL=LLMINfLLMAX 

INDEX  = INDEX  - NMKM1 
76  SCRTCHf INDEX)  = SCRTCH ( LL  4N") 

CALL  PSPLVN  < SCRTCH  fKF'1,2fXIJfILEF‘TfSCRTCH(  1 4NV  ) ) 

ID  * MAXO<  0 f ILEFT-KF  KTl > 

LLMIN  * 1 - MINOIOf ILErT-KRKRl > 

DO  27  LL  =LLM IN » KP 1 
IP  * ID  4 1 

27  SCRTCH < I D4 NA ) « SCRTCH< ID4NA ) - SIGN*SCRTCH(LL4HV) 

28  SIGN  = -SIGN 
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CALI.  PANE  AC  I SCRTCH  ( i 4 NC  > > NMK  > KMi\  .KFKiil  .K.  irLAG) 

fiO  10  (45.44  > , iri  AO 

44  F’RINI  644 

64  4 FORMAT (32H  C IN  SFlOFT  IS  NOT  INVFKTIPLr) 

RE  1 URN 

COMFUlF  SCRTCH  (.4Nr.)  = l>  FROM  SCRTCH  (.4NA>  = - A . 

40  PO  46  I-N.7.-1 

46  SCRTCH(  I-MNA)  = SCRTCH  ( I - H NA ) 4 SCRTCH ( I 4 NA ) 

PO  49  1 1 1 r NMt. 

4V  SCRTCH  ( T 4NP ) - SCFcTCH  ( 1 4NA ) # ( TAU ( 1 4 K ) -T  AU  (I))/F10ATK 
compute  scrt^h  '.4np>  * x . 

CAI  l PANSUl’  < SCRTCH  ( 1 4 NO  . NMK  , NMK  . KFKM  I . K . SCRT  CH  ( 1 4 NP  ) ) 

COMPUTC  SCRTCH  ( . 4 NP ) = CHANGE  IN  XI  . NOPITYf  IF  NECESSARY < TO 
C PREVENT  NEW  XI  FROM  MOVING  MORE  THAN  1/3  OF  THE  WAY  TO  ITS 
C NEIGHPORS.  THEN  APP  TO  XI  TO  OPTAIN  NEW  XI  IN  SCRTCH (. 4NX ) . 
PCI.MAX  * 0. 

SIGN  =-  SIGNST 
PO  53  I-l.NMK 

PEI.  - SION4SCRTCHC I4NP) 

PEI. MAX  = AMAXI  ( PEI. MAX  f APS  ( PEL  > ) 

IF  (PEL  .GT.  0. ) GO  TO  51 

PEL  =■  AMAXI  (PEL.  CSCRTCIK  I- 1 4 NX ) -SCRTCH ( I4NX)  )/3.  > 

GO  TO  5? 

51  PEL  = AMIN1 (PEL. (SCRTCH(I414NX)-SCRTCH(I4NX) >/3.  ) 

5?  SIGN  = -SIGN 

53  SCRTCH! T4NX)  = SCRTCIK  I4NX ) 4 PEL 

CALL  IT  A PAY  IN  CASE  CHANGE  IN  XI  WAS  SMALL  ENOUGH  OR  TOO  MANY 
C STEFS  WERE  TAKEN. 

IF  ( pr  LMAX  .LT.  TOL)  GO  TO  54 

NEWTON  - NEWTON  4 1 

IF  < I4EUT0N  .LE.  NEWTMX ) GO  TO  20 

PRINT  653. NEWTMX 

653  FORMAT  (331.1  NO  EONUERGENCF  IN  SFLOFT  AFTER. 13. 14H  NEWTON  STEFS.) 

54  PO  55  I-l.NMK 

55  T(K4I>  = SCR  TCH (14  NX ) 

56  PC  5 7 1 =■  1 . K 

T ( I > = TALK  1 ) 

57  T(NII)  «>  TAII(N) 

RETURN 

999  IFl  AG  * 2 

RETURN 

END 
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The  subroutine  SPLOPT  has  input  TAU(i)  = x^,  i = l,...,n, 

assumed  to  be  nondecreasing  and  to  satisfy  x^  < x , all  i»  the 

integer  N = n and  the  desired  order  k in  K.  The  routine  needs  a 
work  array  SCRTCH,  of  size  > (n  - k)  (2k  + 3)  + 5k  + 3;  n(2k  + 3) 
is  more  than  enough.  The  routine  has  output  T(i)  = t^, 

i = l,...,n  + k,  the  knot  sequence  for  the  optimal  recovery  scheme, 
in  case  IFLAG  = 1.  For  IF1.AG  = 2,  something  went  wrong. 

The  routine  uses  the  subroutine  BSPLVN  for  the  evaluation  of 
all  B-splines  of  a given  order  on  a given  knot  sequence  which  do 
not  vanish  at  a given  point.  This  routine,  and  others  for  dealing 
computationally  with  splines  and  B-splines,  can  be  found  in  [1]. 

For  completeness,  we  also  list  here  the  subroutines  BANFAC  and 
BANSUB,  used  in  SPLOPT  for  the  solution  of  the  banded  system  (20). 

SUBROUTINE  PANFAC  ( A.  NROW,  N,  NPIAG.  NIPPLE,  IFLAG  ) 

[i  I HT  NS  I ON  A (NROW, NPIAG) 

ITLAG  = 1 

ILO  - NIPPLE  - 1 

IF  (ILO)  * - 999 f 1 0 » 1 9 

10  P0  11  1 = 1 »N 

IF  (ACIr 1))  1 1 » 999 ,11 

11  CUNTINUE 

RETURN 

19  IHI  = NPIAG  - NIPPLE 

IF  (IHI)  999 » 20 , 29 

20  on  25  1 = 1.  N 

IF  ( A ( I » MI  POLE ) ) 21.999.21 

21  JMAX  = NIN0( ILO.N-I ) 

IF  ( JNAX ) 25.20.22 

22  P0  23  J= 1 » JNAX 

23  A( I+J.NIPPLF-J)  = A(I+J»NIPPLE-J)/A< I. NIPPLE) 

25  CONTINUE 

RETURN 

29  DO  50  1 = 1 .N 

PIAG  = A(I, NIPPLE) 

IF  (PIAG)  31,999,31 

31  JNAX  = NTN0( IL0.N-I ) 

IF  (JNAX)  50.50,32 

32  KNAX  = NIN0( IHI »N-I ) 

PO  33  J= 1 . JNAX 

IPJ  = 14  J 
HMJ  = NIPPLE-J 

A(  If'J.NNJ)  = A ( IP J» MN J ) / PIAG 
DO  33  K-l »RNAX 

33  A(  IPJ.MNJ  4 N ) = A ( IF' J , HH  JfN  ) - A ( IP  J , NN  J ) *A  ( I , HIPPLE+N ) 

50  CONTINUE 

RETURN 

999  IFLAG  = 2 

RETURN 

FNP- 
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BANFAC  factors  an  N * N band  matrix  C whose  NDIAG  bands  are  con- 
tained in  the  columns  of  the  NROW  * NDIAG  array  A,  with  the  MIDDLE 
column  containing  the  main  diagonal  of  C.  It  uses  Gauss  elimina- 
tion without  pivoting  and  stores  the  factors  in  A. 


SUBROUTINE  PANSUB  ( A . NROU.  N.  NDIAG.  MIDDLE.  B ) 
DIMENSION  A NROW . NDIAG  ) . D ( N ) 

IF  (N  .Ed.  1)  GO  TO  21 

ILO  = MIDDLE  - 1 

IF  < I L 0 ) 21.21.11 

11  DO  19  1=2 . N 

JMAX  = MINCK  1-1 . ILO) 

DO  19  J=1 » JMAX 

19  P(I)  = B < I ) - B(I-J)*A(I » MIDDLE- J) 


21  I = N 

IH"  = NPIAG-MIDDLE 
DO  29  11=1. N 

JMAX  = MINO (N-I.IHI) 

IF  (JMAX)  28.28.22 

22  DO  25  J= 1 . JMAX 

25  PCI)  = PCI)  - P<I+J)*A<I.MIDDLE+J) 

28  P(I)  = P(I)/A(I. MIDDLE) 

29  1=1-1 
END 


BANSUB  then  uses  the  factorization  of  C into  the  product  of  a 
lower  and  an  upper  triangular  matrix  computed  in  BANFAC  to  solve 
the  equation  Cx  = b for  given  b (input  in  B)  by  forward  and  back 
substitution.  The  solution  x is  contained  in  B,  on  output. 


6.  CONSTRUCTION  OF  THE  OPTIMAL  INTERPOLANT 


With  the  break  points  < ...  < Er  ^ for  the  optimal  inter- 

polant  ?f  determined  from  x in  SPLOPT,  it  remains  to  compute  Sf. 
This  we  propose  to  do  by  determining  its  B-spline  coefficients. 

n+k 

SPLOPT  has  produced  the  knot  sequence  t = (t^)^  , with 


V Vi 


V 1 


1, 


,n  - k, 
. t 


n+1  • ’ ’ "n+k 

Let  (N^)"  be  the  corresponding  sequence  of  normalized  B-splines 
of  order  k,  i.e. , 

Nt(x)  Ni,k,t(x)  (ti+k"ti^ti,‘,,,ti+k^*~x)+  » a11  1- 
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Then,  according  to  Curry  & Schoenberg  [7],  every  piecewise  poly- 
nomial function  of  order  k on  [a,b]  :=  [t,,t  ],  with  k - 2 

1 n 

continuous  derivatives  and  break  points  £^,...,5  i.e.,  every 

spline  of  order  k on  [a,b]  with  knot  sequence  t,  has  a unique 

representation  as  a linear  combination  of  the  n functions 

N N . Therefore 

i n 


Sf 


n 

= l 

i=l 


aiNi 


with  a = (a^)^  the  solution  of  the  linear  system 
n 

(23)  l N (t  )a.  = f(x  ),  i = l,...,n  . 
j-1  J J 


In  case  t is  st 
Lemma  1 implies 


rictly  increasing,  - the  only  case  considered  here, 
that  t^  < < t^+^,  all  i,  which,  together  with 


the  fact  that  N.  vanishes  outside  the  interval  [t  ,t  , ],  all  j, 
J . J 3+'t 


shows  that  the 
the  coefficient 
(see  [5])  solve 
within  the  2k  - 
entries  of  the 
the  linear  syst 


coefficient  matrix  of  (23)  is  2k  - 1 banded.  Since 
matrix  is  also  totally  positive,  we  can  therefore 
(23)  by  Gauss  elimination  without  pivoting  and 
1 bands  required  for  the  storage  of  the  nonzero 
matrix.  ’The  following  subroutine  SPLINT  generates 


em  (23),  given  on  input  the  arrays  TAU(i)  = t , 
FTAU(i)  = f^),  i = 1 N,  T(i)  = t^  i = 1,...,N  + K andthe 


numbers  N and  K.  The  system  is  then  solved,  using  BANFAC  and 
BANSUB  given  in  Section  5,  and  using  a working  array  Q,  of  size 
N(2K  - 1).  The  output  consists  of  the  B-coefficients 
ai  - BCOEF(i),  i = 1 N,  in  case  I FLAG  = 1.  If  I FLAG  = 2,  then 

the  linear  system  (23)  was  not  invertible. 
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SUPR'OUTI  HI  GPLJNT  < TAU.  F TAll  > T»  N>  K,  Cl.  PCOEF.  IFLAG  ) 
c SFL 1 Ml  PRODUCTS  THE  D SPLINE  COrFF.S  PCOEF  OF  1 HE  SPLINE  OF  ORDER 

C K WITH  KNOIS  T (I),  I - 1 ....  . N -I  K . WHICH  TAKES  (IN  THE  VALUE 

C FTAU  (I)  AT  TAU  <I).  1=1....*  N . 

C TAU  IS  ASSUMFH  TO  I<E  STRICTLY  INCREASING. 

C SEE  TEXT  FOR  DESCRIPTION  OF  VARIABLES  ANLi  METHOD. 

C DIMENSION  T (NIK ) »G<N,2*K-1 > 

DIMENSION  TAU(N),  FTAUCN),  TC1>»  QCN.l).  DCOEFCN) 

NF'l  = N + 1 
KF'KMl  = 2*K  - 1 
ILF  FT  = K 
DO  30  I = 1 » N 

TAUI  = TALK  I) 

I LP 1 MX  - M 1 NO ( I 4K  » NF'l  ) 

DO  13  J=1  .KF’KMl 
13  O(I.J)  - 0. 


ILEFT  = MAXO < ILEFT , I ) 

IF  (TAUI  .LT.  T ( ILEFT ) ) 

GO 

TO 

998 

15 

IF  (TAUI  .LT.  T ( ILEFT  + 1 ) ) 

GO 

TO 

16 

ILEFT  * ILEFT  4 1 
IF  (ILEFT  .LT.  ILF1MX) 

GO 

TO 

15 

ILEFT  = ILEFT  - 1 
IF  (TAUI  .GT.  T ( ILEFT  + 1 ) ) 

GO 

TO 

998 

16  CALL  PSF'LVN  ( T.  K»  1.  TAUI.  ILEFT , PCOEF  ) 

C NOTE  THAT  PCOEF  IS  USED  HERE  FOR  TEMP . STORAGE . 

L = ILEFT  - I 
DO  30  J=1.N 
L = L + l 

30  0(1, L)  = PCOEF(J) 

NF'2MN  = N42-K 

CALL  PANFAC  ( 0.  N.  N»  KF'KMl . K.  IFLAG  ) 

GO  TO  (40.999).  IFLAG 

40  DO  41  1 = 1,  N 

41  PCOEF ( I ) = FTAU ( I ) 

CALL  PANSUP  C 0,  N,  N.  KF'KMl,  K,  PCOEF  ) 

RETURN 

998  IFLAG  = 2 

999  PRINT  699 

699  FORMAT  < 4 1 H LINEAR  SYSTEM  IN  SPLINT  NOT  INVERTIPLE) 

RETURN 

END 


Note  that  [1]  contains  programs  w Hch  might  facilitate  sub- 
sequent use  of  the  optimal  interpolant  determined  in  this  way 
via  SPLOPT  and  SPLINT. 

Finally,  the  linear  system  (23)  can  be  generated  in  0(nk2) 
operations  and,  because  of  the  band  structure,  can  be  solved  in 
0(nk)  operations.  The  linear  system  (20),  to  be  generated  and 
solved  at  each  Newton  step  for  finding  £,  is  of  similar  nature 
(with  n - k rather  than  n equations  and  a coefficient  matrix  which 
is  the  transposed  of  the  kind  of  matrix  appearing  in  (23))  hence 
requires  a similar  effort  for  its  construction  and  solution. 
Therefore,  if  it  takes  indeed  only  three  to  four  Newton  iterations 
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to  find  C to  sufficient  accuracy,  then  It  takes  only  four  to  five 
times  as  much  work  to  construct  the  optimal  interpolant  rather 
than  any  spline  interpolant  to  the  same  data.  Also,  the  total 
effort  is  only  0(nk2)  which,  for  large  n,  compares  very  favorably 
with  such  interpolation  schemes  as  polynomial  interpolation  which 
takes  O(n^)  operations,  or  more  general  schemes  which  take  as 
much  as  O(n^)  operations. 
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